Numerical study of the disorder-driven roughening transition in an 
elastic manifold in a periodic potential 
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We investigate the roughening phase transition of a (3 + l)-dimensional elastic manifold driven 
by the completion between a periodic pinning potential and a randomly distributed impurities. The 
elastic manifold is modeled by a solid-on-solid type interface model, and universal features of the 
transition from a flat phase (for strong periodic potential) to a rough phase (for strong disorder) 
are studied at zero temperature using a combinatorial optimization algorithm technique. We find a 
continuous transition with a set of numerically estimated critical exponents that we compare with 
analytic results and those for a periodic elastic medium. 



PACS numbers: 64.60.Cn, 68.35.Ct, 75.10.Nr, 02.60.Pn 



Extended objects embedded in a higher-dimensional 
space, like polymers, magnetic flux lines, surfaces, inter- 
faces or domain walls, are on large length scales com- 
monly described by elastic models, the so-called elastic 
manifolds (EM) Q| . If quenched disorder in form of im- 
purities or other randomly distributed pinning centers is 
present, the EM will be, even at low temperatures or in 
the absence of any thermal fluctuations, in a rough state. 
However, a periodic array of pinning potentials, like a 
background lattice potential, increases the tendency of 
the EM to minimize their elastic energy, i.e. to stay flat 
(3). Both mechanism compete with each other and 
by varying their relative strength a roughening transi- 
tion might emerge. The numerical investigation of such 
a scenario is the purpose of this paper. 

Consider a d-dimcnsional EM in a (d+ l)-dimensional 
medium with quenched impurities distributed randomly. 
Fluctuations of the shape of the EM are then described 
by a scalar displacement field </>(r) denoting a deviation 
from a flat reference state in the (d+l)th direction at each 
r e R d ; (</>(r), r) refers to the (d+ l)-dimensional coordi- 
nate of an EM segment. It is known that quenched disor- 
der, no matter how weak, destabilizes the flat phase for 
d < 4 H . The emerging disordered rough phase is charac- 
terized by a divergent displacement correlation function 
B(r) = ([<p(r) — </>(0)] 2 ) with the distance r = |r|, whose 
scaling property is universal. 

An example of the ID EM is a directed polymer or 
a magnetic flux line in a disordered 2D plane. Using a 
mapping to the Kardar-Parisi-Zhang equation for surface 
growth, it is shown analytically that the displacement 
correlation diverges algebraically as 



B(r) 



(1) 



with the exactly-known roughness exponent £ = | [Q. 
For higher-dimensional EMs, analytic studies using a 
functional renormalization group (FRG) method predict 
that the rough phase is governed by a zero-temperature 
fixed point which is also characterized by a power-law 
divergence of B(r) as in Eq. (|l|) ||, || and the rough- 
ness exponent is found to be £ = 0.2083e + C(e 2 ) § 
or |e + 0{e 2 ) || up to first order in e — 4 — d. These 



values are in good agreement with numerical estimates 
obtained from exact ground-state calculations for (2 + 1) 
and (3 + 1) D solid-on-solid type interface models (j], §]. 

In addition to random impurities, also the structure of 
the embedding medium can affect the large scale proper- 
ties of the EM. In particular, if the medium has a crys- 
talline structure, the EM is pinned by the disorder po- 
tential and by the periodic crystal potential, both effects 
competing with each other — the latter favoring a flat 
state while the former a rough one. Hence, the EM may 
undergo a phase transition at a critical disorder-potential 
strength. Following a qualitative perturbative scaling ar- 
gument ||, one can see that there exists the disorder- 
driven roughening transition at nonzero disorder strength 
for 2 < d < 4. 

In 2D, the disorder potential is argued to dominate 
over the periodic potential marginally ||, ||. Conse- 
quently the 2D EM is believed to be rough at any dis- 
order strength asymptotically beyond a certain length 
scale which diverges exponentially with the inverse of a 
disorder strength. Some numerical studies support such 
claim |)J. On the other hand, it has been reported that 
an Ising domain wall in a (2+l)D lattice with bond dilu- 
tion displays a roughening transition at a non-zero dilu- 
tion probability B . It suggests that the type of disorder 
might be important in the marginal 2D case (see also 
discussions in Ref. M). 

In 3D, the existence of a roughening transition was 
shown in the studies using a Gaussian variational (GV) 
method § and a FRG method g, In the GV 

study, the free energy was calculated by approximating 
the Hamiltonian of the EM with a Gaussian. It leads to 
a conclusion that the transition is of first order. On the 
other hand, the FRG study with a perturbative expan- 
sion in the periodic potential strength and in e = 4 — d 
yields that the transition is continuous with a correlation 
length exponent v = l/(2 v / e). To clarify this issue we 
performed in this work a numerical study of the disorder- 
driven roughening transition of the 3D elastic manifold 
in a crystal potential with quenched random impurities. 
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The EM is described by the Hamiltonian 



H= d d 



: |V0(r)| 2 +F P (,/)(r)) + ^( ( /»(r),r) , (2) 



where the first term represents a surface tension, Vp(<j)) = 
— Vcoscj) the periodic lattice potential, and Vr{4>, r) the 
quenched random potential. Here <fi is measured in units 
of ao/(27r) with ao the lattice constant of the crystal. For 
uncorrelated distribution of impurities, Vr(0, r) can be 
taken as a random variable with mean zero and variance 
given by 



V R ((/>, r)V R (<P', r') = D R(<j> - 0') S(r - r') (3) 

with a parameter D for the disorder strength. Uncor- 
related distribution of impurities implies that the disor- 
der correlation function in the (d + l)th direction is also 
short-ranged, i.e., R(<fi) = S(cj)). 

It is interesting to note that the Hamiltonian in Eq. (|^) 
can also describe so-called periodic elastic media (PEM) 
in a crystal with quenched disorder ||. A system of 
strongly interacting particles or other objects, like mag- 
netic flux lines in a type-II superconductor or a charge 
density in a solid, will order at low temperatures into a 
regular arrangement, namely, the flux line lattice (FLL) 
or the charge density wave (CDW), respectively. Fluc- 
tuations either induced by thermal noise or by disorder 
induce deviations of the individual particles from their 
equilibrium positions. As long as these fluctuations are 
not too strong, an expansion of the interaction energy 
around these equilibrium configuration might be appro- 
priate. An expansion up to second order leads to the 
surface-tension- like term as in Eq. ^ . In contrast to the 
EM, the PEM have their own periodicity A. It implies 
that the disorder potential Vr(4>,t) should be a periodic 
function in <fi with the periodicity p = ^ (commensu- 
rability parameter), even though the impurities are dis- 
tributed randomly. Hence, the disorder correlation func- 
tion in Eq. (||) should be periodic: R(cj> + 2np) = R{4>). 
In 3D, as a result of the periodicity, the displacement 
correlation function for tp = <f>/p diverges logarithmically 
as ([(fi(r) — ip(0)] 2 ) ~ 2A\nr with a universal coefficient 
A ~ 1.0 in the rough phase Jll], [l^, The periodic 
elastic media also display a disorder-driven roughening 
transition as a result of competition between the periodic 
potential and the random potential. However, there is a 
slight controversy regarding the nature of the transition 
since an analytic FRG study and a zero-temperature 
numerical study yield results that are not fully com- 
patible. 

We introduce a discrete solid-on-solid (SOS) type in- 
terface model for the elastic manifold whose continuum 
Hamiltonian is given in Eq. @ . Locally the EM remains 
fiat in one of periodic potential minima at 4> — 2ixh with 
integer h. Due to fluctuations, some regions might shift 
to a different minimum with another value of h to create a 
step (or domain wall) separating domains. To minimize 
the cost of the elastic and periodic potential energy in 



Eq. (||), the domain- wall width must be finite, say £ Q ||. 
Therefore, if one neglects fluctuations in length scales less 
than £ , the continuous displacement field cf>(r) can be re- 
placed by the integer height variable {/i x } representing 
a (3 + 1) D SOS interface on a simple cubic lattice with 
sites x £ {1, . . . , L} 3 . The lattice constant is of order 
and set to unity. The energy of the interface is given by 
the Hamiltonian 



H = 



J(h x ,x);(h y ,y)\hx- hy\ -^Vr(/i x ,x) , (4) 

where the first sum is over nearest neighbor site pairs. 
After the coarse graining, the step energy J > as well 
as the random pinning potential energy Vr becomes a 
quenched random variable distributed independently and 
randomly. Note the PEM has the same Hamiltonian as in 
Eq. (0) with random but periodic J and Vr in h with peri- 



odicity p 13 . In this sense, the elastic manifold emerges 
as in the limit p — > oo of the periodic elastic medium. 

Here, we are interested only in the ground state prop- 
erty. Since the prevalent RG-picture suggests that the 
roughening transition is described by a zero temperature 
fixed point || 10, |Il| one expects the critical exponents 
we find to be valid for the finite temperature roughening 
transition as well. To find the ground state, one maps the 
3D SOS model onto a ferromagnetic random bond Ising 
model in (3 + 1)D hypercubic lattice with anti-periodic 
boundary conditions in the extra dimension Q (for the 
3 space direction one uses periodic boundary conditions 
instead). The anti-periodic boundary conditions force a 
domain wall into the ground state configuration of the 
(3+l)D fcrromagnet. Note that bubbles are not present 
in the ground state. A domain wall may contain an over- 
hang which is unphysical in the interface interpretation. 
Fortunately, one can forbid overhangs in the Ising model 
representation using a technique described in Ref. [Q. 
If the longitudinal and transversal bond strengths arc 
assigned with J/2 and Vr/2 occurring in Eq. ([|), re- 
spectively, this domain wall of the ferromagnet becomes 
equivalent to the ground state configuration of ([|) for 
the interface with the same energy. The domain wall 
with the lowest energy is then determined exactly by us- 
ing a combinatorial optimization algorithm, a so-called 
max-flow/min-cost algorithm. This combinatorial opti- 
mization technique is nowadays standard in the study of 
disordered systems and we refer readers to Ref. [[l5| for 
a detailed review. 

We performed the ground state calculation on L 3 x H 
hypercubic lattices for L < 32. H, the size in the extra 
direction, is taken to be larger than the interface width. 
Several distributions for J and Vr were studied for the 
critical behavior of a disordered system may depend on 
the choice of the disorder distribution as in the random 
field Ising system Jl6| |. However, our main numerical re- 
sults do not depend on the specific choice of the distri- 
bution. So we only present the results for an exponential 
distribution for J > 0, P(J) 
distribution for < Vr < V ma . 
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FIG. 1: (a) W vs. L in a log-log plot at A = 8 and 10 
(strong disorder), (b) £ vs. 1/L for the same values of A 
as in (a). The solid lines are least square fits to the form 
C(L) = ( + aL' b . We obtain £ = 0.21 ±0.01 and the resulting 
value of 6 ~ 1.5 > 1 indicates that the extrapolation is stable 
against statistical uncertainties. 



controlled with the parameter A = V max / Jo- Other dis- 
tributions studied include (bimodal,bimodal) and (uni- 
form, uniform) distributions for (J, Vr), and gave identi- 
cal estimates for the critical exponents. 

The state of the i nterface is char acterized by the width 
W defined as W 2 = (h%) — (/i x )o , where (■ • -) denotes 
the spatial average in the ground state and ~ the disor- 
der average. W 2 is proportional to the spatial integral of 
the displacement correlation function. We also measure 
a magnetization-like quantity m = \ (e l7rh *) \. It is analo- 
gous to the magnetization used as an order parameter for 
the roughening transition of the PEM [jl3| . One expects 
that to is non-zero in the flat phase and vanishes in the 
rough phase. So it can be used as an order parameter for 
the roughening transition. 

We first examine the power-law scaling behavior of 
the width, W ~ Z>, in the rough phase at large dis- 
order strength A = 8 and 10. Figure |l| (a) shows the 
width average over 1000 ~ 5000 disorder realizations 
for L — 4 ~ 32. A noticeable curvature in the log- 
log plot indicates that corrections to scaling are still 
rather strong. Nevertheless, we can estimate the rough- 
ness exponent by extrapolating an effective exponent 
C(L) = \og[W(2L)/W(L)]/\og2 to the limit L -» oo by 
fitting it to a form ((L) = £ + aL~ fc , see Fig. (b). We 
obtain that 



C = 0.21 ±0.01 



(5) 



This value is consistent with the previous analytical and 
numerical results ||, Q, || . 

As the disorder strength decreases, the width also de- 
creases and the interface eventually becomes flat below a 
threshold (see Fig. |J). The order parameter also shows 
an indication of a phase transition around A ~ 4.0. 
Apparently the order parameter decreases continuously. 
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FIG. 2: The interface width (a) and the order parameter (b) 
as functions of A at L — 4, 8, 16, and 32. They are obtained 
from the disorder average over 3000 ~ 20000 samples. 
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FIG. 3: The effective exponent [/3/u](L) for different values 
of A as a function of 1/L. The broken line is a guide for 
the eyes that separates the curves with a downward bending 
(A < A c ) from those with an upward bending (A > A c ). 



Therefore we perform a scaling analysis assuming that 
the phase transition is a continuous one. The critical 
point A c can be determined from the finite-size-scaling 
property of the order parameter: 



m(L, e) 



-fi/ v 



J={eL x l v ) 



(6) 



where e = A — A c , and (3 (is) is the order parame- 
ter (correlation length) exponent. The scaling function 
T{x) has a limiting behavior T{x — > 0) = const, so 
that the order parameter decays algebraically with L 
as to ~ L~PI V at the critical point. It also behaves as 
T(x — * — oo) ~ \x\P so that to ~ \z\& for A < A c in the 
infinite system size limit. Consider the effective expo- 
nent \f)/v](L) = — log(m(2L)/m(L))/log2. It converges 
to the value of @/u at the critical point and deviates 
from it otherwise as L increases. We estimate the criti- 
cal threshold as the optimal value of A at which the ef- 
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FIG. 4: Scaling plot of mi"/" vs. \e\L 1/v with e = A — 4.25, 
= 0.07, and v = 1.4. 



30 p 

- 20- 
g 

^ 10- 



0^ 



0.2 



6;,, 



0.4 0.6 

m 



0.8 




A = 4.0 
A = 4.2 
A = 4.4 
A = 4.6 
A = 4.8 



A = 5.0 
A = 5.2 
A = 5.4 
A = 5.6 
A = 5.8 



FIG. 5: Histograms for the probability distribution P(m) of 
the magnetization at L — 32. 



fective exponent approaches a nontrivial value. The plot 
for this effective exponent is shown in Fig. ||. One can 
see that there is a downward and upward curvature for 
A < 4.20 and A > 4.30, respectively. From this behavior 
we estimate that A c = 4.25 ± 0.05 and 

^=0.07 ±0.03. (7) 

Note that the effective exponent varies with L even at the 
estimated critical point, which implies that corrections to 
scaling are not negligible for system sizes up to L = 32. 
For that reason our numerical results for A c and (3/v have 
rather large error bars, and one may need larger system 
sizes for better precision. The exponents (3 and v could 
also be obtained from the scaling analysis using Eq. (||). 
We fix the values of A c and (5/v to the values obtained 
before and vary v to have an optimal data collapse. We 
obtain 

i/ = 1.4±0.2 (8) 

and the corresponding scaling plot is shown in Fig. ||. 

The order-parameter scaling property shows that the 
roughening phase transition is a continuous transition, 
though the exponent (3 ~ 0.1 is very small, as opposed 
to the results of the GV study j!| predicting a first order 
transition. The transition nature becomes more trans- 
parent by looking at the probability distribution P(m) 
of the magnetization near the critical point. We mea- 
sured the distribution from a histogram of the magneti- 
zation of 3000 samples with L = 32, which is shown in 
Fig. [j| We do not observe a double-peak structure in 
P(m), which would appear for a first order transition, at 
any values of A. Instead, there is a single peak which 
shifts continuously toward zero as A increases. We did 
not observe any double-peak structure in the distribu- 
tion of the width, either. Therefore we conclude that the 
roughening transition is the continuous transition. 



We note that this behavior is reminiscent of the three- 
dimensional random field Ising-model (RFIM) [jl7| , where 
also a very small order is found. This causes an extremely 
weak systems size dependence and the transition appears 
to be discontinuous in the order parameter although the 
system is indeed critical and the correlation length di- 
verges. However, as has been pointed out in Jl8| the 
magnetization of an individual sample shows close to the 
transition discontinuous jumps when (for a fixed field dis- 
tribution) the coupling strength is varied (however, the 
sample-averaged magnetization is smooth). Since in the 
3d RFIM the maximum jump does not vanish even in 
the infinite size limit some objection against the continu- 
ity of the phase transition has been raised . Since we 
find also a small order parameter exponent - and because 
of some model specific similarities between the 3d RFIM 
and the system we study here — we now want to check 
the jump size statistics, too. 

For a given realization of J and Vr in Eq. (|J) , we mea- 
sure the order parameter and its jumps when changing 
Vr with a global factor. The interface may undergo two 
types of intermittence: The average height (h) may jump 
with a vanishing overlap of the interface configurations 
before and after the jump (meaning that the whole man- 
ifold is in a new position, uncorrelated with the previous 
one) . Or a large domain type excitation may appear with 
only a small change in (h) . The intermittent behavior of 
(1+1) and (2+l)D interfaces were studied in great detail 
in Ref . || . We want to study how the interface at a given 
average height roughens, therefore we only take into ac- 
count the domain type excitations with the change in (h) 
less than one half when measuring the jump of the order 
parameter. 

We calculate the order parameter in the interval 3 < 
A < 7 with spacing 0.02. Figure ^| shows the probability 
distribution of maximum values of the jump 5m in 1000 ~ 
5000 samples of sizes L = 4, 6, 8, 12, 16. The inset shows 
the averaged value. It is non-zero, but becomes smaller 
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FIG. 6: Order parameter jump distribution. Inset shows the 
plot of the averaged value vs. 1/L in the log-log scale. The 
solid line has a slope of 0.16. 



FIG. 7: (a) Effective exponent for the power-law behavior 
of W ~ L . (b) Prefactor in the logarithmic scaling of W 2 ~ 
AlogL. 



as the system size increases except for the case of L = 4. 
The decay is very slow, nevertheless we could fit the data 
for L > 6 to a form Sm ~ L~ 016 as can be seen in 
the inset. It suggests that the order parameter jump 
vanishes in the infinite size limit, and the smallness of 
the estimated exponent 0.16 is compatible with our small 
estimate for the order parameter (3. 

We also studied the scaling of the width at the critical 
point. One might expect that the width at the critical 
point scales as W ~ with a new roughness exponent 
£' different from £ ~ 0.21 for the rough interface. We 
plot the effective exponent C = log[W(2£)/W(£)]/log2 
near the critical point in Fig. [?] (a). One observes that 
this effective exponent does not extrapolate to a non- 
zero value (for comparison c.f. Fig. 0(b)). Instead, it 
decreases rapidly as L grows, which suggests rather a 
logarithmic scaling, W 2 ~ AlogL, at the critical point 
as in the case of the periodic elastic medium |r|. To 
investigate such a possibility, we measure the prefactor 
A = (W 2 (2L) - VF 2 (L))/log2 for this logarithmic scal- 
ing. They are also plotted in Fig. fj] (b), which shows 
a clear threshold behavior: It decreases (increases) for 
A < 4.25 (> 4.25) and remains constant (A ~ 0.03) 
at A = 4.25, which was estimated as the critical point 
from the order parameter scaling analysis. These facts 
are consistently suggesting that the interface width scales 
logarithmically at the critical point: 



W 2 ~ 0.03 log L 



(9) 



The result of the logarithmic scaling of the width at the 
critical point agrees well with that of the FRG study jio) . 

Our numerical results pose an interesting question: In 
recent work fl3|| we found that the disorder-driven rough- 
ening transition of the PEM appears to be independent 
of the commensurability parameter p jl3| . As mentioned 
before, the elastic manifold can be seen as the p — > oo 



limit of the PEM. Based on this observation one might 
speculate that the roughening transitions of both sys- 
tems belong to the same universality class. Although the 
p — > oo-limit of the PEM could also belong to a different 
universality class the numerical results are compatible 
with the roughening transition in the EM and the PEM 
belonging to the same universality class:: For the EM we 
obtain j3/u ~ 0.07, v ~ 1.4, and A ~ 0.03 (W 2 ~ AlogL 
at the critical point), while for the PEM [13j we got 
/3/v ~ 0.05, v ~ 1.3, and A ~ 0.018. Although these val- 
ues are very close to each other, a final conclusion could 
not be drawn yet due to the rather strong correction to 
finite size scaling observed in our study of the EM. So we 
have to leave this issue as an open question. 

In summary, we have studied the (3 + 1)D elastic man- 
ifold in a crystal with quenched random impurities. We 
have investigated numerically the disorder-driven rough- 
ening transition at zero temperature using an exact com- 
binatorial optimization algorithm technique. The tran- 
sition turns out to be continuous with the critical ex- 
ponents (3/v ~ 0.07 and v ~ 1.4. For a given disor- 
der potential configuration, the order parameter shows 
a discrete jump for finite size systems when varying its 
strength. However, the jump vanishes in the infinite sys- 
tem size limit. We also found that W 2 scales logarithmi- 
cally with the system size L at the critical point, in con- 
trast to the power-law scaling in the rough phase. Our 
results do not agree with the scenario proposed in Ref. 
that the roughening transition should be first order. In- 
stead, they are in a qualitative agreement with those of 
FRG e-expansion study in Ref. flTo{ , which predicts a 
continuous roughening transition and a logarithmic di- 
vergence of W 2 at the critical point. However, there is a 
significant discrepancy between the values of the critical 
exponents obtained numerically and analytically, which 
was also observed in the study of the PEM jl3) . 
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